home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / testmatrix / pnorm.r < prev    next >
Encoding:
Text File  |  1994-12-27  |  3.2 KB  |  137 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Estimate of matrix p-norm (1 <= p <= inf).
  4.  
  5. // Syntax:    pnorm ( A , p , TOL ) 
  6.  
  7. // Description:
  8.  
  9. //    Estimates the Holder p-norm of a matrix A, using the p-norm
  10. //    power method with a specially chosen starting vector. TOL is a
  11. //    relative convergence tolerance (default 1E-4). 
  12.  
  13. //    Returned as list elements are the norm estimate EST (which is
  14. //    a lower bound for the exact p-norm), the corresponding
  15. //    approximate maximizing vector x, and the number of power
  16. //    method iterations k. A nonzero fourth argument causes trace
  17. //    output to the screen. 
  18.  
  19. //    If A is a vector, this routine simply returns NORM(A, p). 
  20.  
  21. //      See also NORM, NORMEST.
  22.  
  23. //      Note: The estimate is exact for p = 1, but is not always exact
  24. //    for p = 2 or p = inf.  Code could be added to treat p = 2 and
  25. //    p = inf separately.
  26.  
  27. //      Calls DUAL and SEQA.
  28.  
  29. //       Reference:
  30. //       N.J. Higham, Estimating the matrix p-norm,
  31. //       Numer. Math., 62 (1992), pp. 539-555.
  32.  
  33. //    This file is a translation of pnorm.m from version 2.0 of 
  34. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  35. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  36.  
  37. // Dependencies
  38.    require dual seqa
  39.  
  40. //-------------------------------------------------------------------//
  41.  
  42. pnorm = function ( A , p , tol , noprint )
  43. {
  44.   local ( A , p , tol , noprint )
  45.   global (pi)
  46.  
  47.   if (!exist (p)) { 
  48.     error ("pnorm: Must specify norm via second parameter.");
  49.   }
  50.  
  51.   m = A.nr; n = A.nc;
  52.   if (min(m,n) == 1) {
  53.    return norm (A, p);        // Return vector p-norm 
  54.   }
  55.  
  56.   if (!exist(noprint)) { noprint = 0; }
  57.   if (!exist(tol)) { tol = 1e-4; }
  58.  
  59.   // Stage I.  Use Algorithm OSE to get starting vector x for power method.
  60.   // Form y = B*x, at each stage choosing x(k) = c and scaling previous
  61.   // x(k+1:n) by s, where norm([c s],p)=1.
  62.  
  63.   sm = 9;    // Number of samples.
  64.   y = zeros(m,1); 
  65.   x = zeros(n,1);
  66.  
  67.   for (k in 1:n)
  68.   {
  69.     if (k == 1)
  70.     {
  71.       c = 1; s = 0;
  72.     else
  73.       W = [A[;k], y];
  74.  
  75.       if (p == 2)    // Special case.  Solve exactly for 2-norm.   
  76.       {
  77.         Wsvd= svd(W);    // [U,S,V]
  78.         c = Wsvd.vt'[1;1]; s = Wsvd.vt'[2;1];
  79.       else
  80.         fopt = 0;
  81.         for (th in seqa(0,pi,sm))
  82.         {
  83.           c1 = cos(th); s1 = sin(th);
  84.           nrm = norm([c1, s1], p);
  85.           c1 = c1/nrm; s1 = s1/nrm;    // [c1, s1] has unit p-norm.
  86.           f = norm( W*[c1, s1]', p );
  87.           if (f > fopt)
  88.           {
  89.             fopt = f;
  90.             c = c1; s = s1;
  91.           }
  92.         }
  93.       }
  94.     }
  95.  
  96.     x[k] = c;
  97.     y = x[k]*A[;k] + s*y;
  98.     if (k > 1) { x[1:k-1] = s*x[1:k-1]; }
  99.  
  100.   }
  101.  
  102.   est = norm(y, p);
  103.   if (noprint) { printf("Alg OSE: %9.4e\n", est); }
  104.  
  105.   // Stage II.  Apply Algorithm PM (the power method).
  106.  
  107.   q = dual(p);
  108.   k = 1;
  109.  
  110.   while (1)
  111.   {
  112.     y = A*x;
  113.     est_old = est;
  114.     est = norm(y, p);
  115.  
  116.     z = A' * dual(y,p);
  117.  
  118.     if (noprint)
  119.     {
  120.       printf("%2.0f: norm(y) = %9.4e,  norm(z) = %9.4e", ...
  121.                  k, norm(y,p), norm(z,q))
  122.       printf("  rel_incr(est) = %9.4e\n", (est-est_old)/est);
  123.     }
  124.  
  125.     if (( norm(z,q) <= z'*x || abs(est-est_old)/est <= tol ) && k > 1)
  126.     {
  127.        return est;
  128.     }
  129.  
  130.     x = dual(z,q);
  131.     k = k + 1;
  132.  
  133.   }
  134.  
  135.   return << est = est; x = x; k = k >>;
  136. };
  137.